One-dimensional continuum electronic structure calculations with the density matrix 

renormalization group 
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We extend the density matrix renormalization group to compute exact ground states of continuum 
many-electron systems in one dimension with long-range interactions. To demonstrate the power 
of this approach, we find the exact ground state of a chain of 100 strongly correlated artificial 
hydrogen atoms. We expect this method to be important for simulating Id cold atom systems and 
for studying density functional theory (DFT) in an exact setting. Not only can we compare DFT 
approximations with exact results, but we can implement the exact density functional. 
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For electronic structure calculations, these are the best 
of times and the worst of times. When correlations are 
weak, density functional theory (DFT) makes it possible 
to tackle extremely realistic Hamiltonians and large sys- 
tem sizes with reasonable accuracy [1]. For strongly cor- 
related systems, there exist powerful and controllable nu- 
merical methods [2] for simulating model lattice Hamil- 
tonians, such as the Hubbard model. However, there are 
few numerical tools that can treat the combination of 
strongly correlated electronic systems and realistic micro- 
scopic Hamiltonians. In the strongly correlated regime, 
DFT approximations are neither systematic nor control- 
lable, often leading to unrestrained parameter multipli- 
cation and empiricism. Model Hamiltonians rely on the 
arbitrary truncation of terms that may be crucial in tip- 
ping the balance between competing phases. Attempts to 
bridge the gap between realistic Hamiltonians and strong 
correlation techniques, such as dynamical mean field the- 
ory coupled to DFT, still contain both arbitrary trunca- 
tions and less than ideal treatment of correlations. 

Here we show that at least in the case of one dimension 
(Id), we can treat both the Hamiltonian and the correla- 
tions essentially exactly, even for very large systems. Our 
approach is based on the density matrix renormalization 
group (DMRG) [3], the most powerful of the strongly 
correlated techniques for Id lattice models. Here we ex- 
tend DMRG to treat continuum electron systems with 
long-range interactions. This new approach retains the 
exponential convergence and near linear scaling with sys- 
tem size of DMRG. As an example, we present a near ex- 
act calculation of a system with 100 strongly interacting 
pseudo-hydrogen atoms. 

A key motivation for this method is to study DFT in 
an exact many-electron setting. For example, we can eas- 
ily compare various DFT approximations with exact re- 
sults. Generically, one dimensional systems have strong 
quantum fluctuations and provide severe challenges for 
DFT. In addition, we can implement the exact energy 
functional of the density, the heart of DFT; previously 
this could only be done for a few electrons. Our DMRG 
approach also provides new ways to characterize realistic 



many-electron quantum systems using quantum informa- 
tion concepts, such as the entanglement entropy across 
bipartitions of the system. Another motivation for study- 
ing this class of systems is that Id continuum Hamilto- 
nians can be realized exactly in cold atom systems [4, 5]. 
Presently, real-space DMRG methods for solid-state 
applications are designed only to work with lattice mod- 
els. Each site of such a model can be thought of as a 
Wannier function centered on an atom. One way of gen- 
eralizing this picture to make the Hamiltonian more real- 
istic is to expand in a set of basis functions; DMRG has 
become a powerful technique for the quantum chemistry 
of small molecules based on this idea [6, 7]. Here we pro- 
ceed in a more flexible direction which does not depend 
on a choice of basis and has optimal scaling of calculation 
time with system size: we represent the continuum in Id 
with a real space grid. The continuum Hamiltonians of 
interest can be written as 

+ / v(x)n(x) + / v cc (x — x 1 ) n(x) n(x') . (1) 

J x Jx,x' 

where ijy a ( x ) creates an electron of spin a at position 
x; v and v co are local and electron-electron potentials, 
respectively; and n is the electron density operator. We 
introduce a grid spacing a, obtaining a discretized Hamil- 
tonian 
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where jl = jjl — I /a 2 



v(J a) and u« = v cc (\i - j\ a). 



The Sij in the last term prevents an unphysical self inter- 
action. For technical reasons we work with open bound- 
aries rather than periodic, and extend the grid well past 
the edge atoms. Finite grid spacing errors can be reduced 
arbitrarily by reducing a; convergence can be accelerated 



by using a higher order discretized derivative. Here, we 
fix a = 0.1 in atomic units. 

A different approach was proposed recently by Ver- 
straete and Cirac [8], who showed how to define and op- 
timize matrix product states directly in the continuum 
limit in the context ol quantum field theories. 

Working efficiently with such a Hamiltonian represents 
an unusual challenge for DMRG. Normally, one is most 
concerned with the number of states per block m needed 
to represent the ground state. The number of sweeps Ng 
needed to converge to the ground state is usually quite 
small (~2— 5) for Id systems. Here, the reverse can hap- 
pen: a small grid spacing 'a' relative to the interatomic 
separation can lead to energy scales that differ by orders 
of magnitude, greatly increasing N$ but not affecting m 
significantly. Fortunately, while convergence with m re- 
flects the entanglement of the system, an inherent prop- 
erty, many approaches can be tried to reduce N$. We 
have found a particularly efficient acceleration approach 
based on a real-space RG procedure which produces a 
supplementary grid with a much coarser spacing and 
lower energy scales such that N$ can be made small; after 
sweeping on this grid we map the wavefunction back onto 
the fine grid for further sweeps. This procedure, which 
will be described in a future publication, introduces no 
additional approximations; it merely reduces the compu- 
tational time to reach a converged ground state. 

Besides the wide range of energy scales, another chal- 
lenge for DMRG is the presence of long-ranged interac- 
tions. The simplest approach, dealing with N g terms 
in a sweep over N g grid points, would scale as Ng. A 
more efficient approach using intermediate operators, as 
used with DMRG for quantum chemistry in a basis set, 
would scale as Ng. We utilize a much more efficient ap- 
proach than either of these by using a representation of 
the Hamiltonian as a matrix product operator (MPO). 
Finite bond dimension MPOs naturally encode exponen- 
tially decaying interactions [9] . Interactions with a power 
law decay can be approximated to high accuracy by fit- 
ting to a sum of AmpO exponentials (usually Ampo < 15 
is sufficient to obtain an accuracy of 10 -5 ) [10]. Thus the 
DMRG calculation time is nominally linear in N g ; specif- 
ically, it is proportional to NsN g NMPO m3 - Deviations 
from a purely linear computational effort with the over- 
all system size can come from a dependence of A5 or 
m with system size; in practice we find that the RG ac- 
celeration procedure keeps A5 small. The behavior of 
in is well-understood for short-ranged model Hamiltoni- 
ans: for non-critical systems, m is independent of sys- 
tem size. For critical systems, with power-law decaying 
correlations, m grows only logarithmically with length. 
Hence, we expect only a slightly worse-than-linear com- 
putational effort with system size in the worst case. 

As a simple setting for exploring strong many-body 
correlation effects, we consider chains of one dimensional 
'soft hydrogen' atoms. Each atom is modeled as a single 




FIG. 1. The exact ground state density of a chain of 100 
widely separated artificial atoms. The total length of the 
system is L = 420 in atomic units (4200 grid sites with a 
spacing of 0.1). The upper panel shows the electron density 
of a central region superimposed with the density at the left 
edge (the dashed blue curve with corresponding x axis at the 
top of the plot). The lower panel shows the electron density 
predicted by DFT within the local density approximation, for 
both restricted and unrestricted spin densities. 



electron in a soft-Coulomb potential well of the form 



"atom (^ = -1/Vx 2 + 



1 



(3) 



where we use atomic units throughout. We 

also include repulsive interactions defined by 
v ee (x — x') — —v & t om (\x — x'\). This potential is a 
standard choice for avoiding complications arising from 
divergent Coulomb interactions and has been used to 
study molecules in intense laser fields [11, 12]. We are 
also fortunate to benefit from the work of Ref. 13, which 
provides a parameterization of the Id local density 
approximation (LDA) for just such an interaction. 

As a demonstration of the power of our approach, we 
display in Fig. 1 the exact ground state density of a chain 
of one hundred artificial atoms with long-range interac- 
tions, which took a few days of computer time on a single 
workstation. Representing the ground state accurately 
required keeping about m = 200 states. The energy er- 
ror in the many-body solution using DMRG is of order 
10~ 4 . The errors arising from the finite grid spacing and 
finite number of exponentials Ampo used to fit v cc are 
larger, of order 10~ 2 , but are well understood and easily 
reduced if necessary. Also shown are both restricted and 
unrestricted DFT calculations using the Id LDA [13]. In 
this system both DFT approaches have substantial er- 
rors. In terms of lattice models, one would represent 
such a system with either a half-filled Hubbard chain or 
with an antiferromagnetic Heisenberg chain. Note that 
both of these models are critical with power law decaying 
spin-spin correlations; a non-critical system would have 



been easier for DMRG. It is not surprising that the local 
DFT methods cannot capture the quasi long-range order 
quantitatively, with the unrestricted LDA giving long- 
range antiferromagnetism (similar to Fig. 2). It is some- 
what surprising that even the total density from DFT 
deviates strongly from the exact DMRG results. The to- 
tal energy from either DFT calculation is off by about 
1%. 

In Fig. 2 we show a system which reveals weaknesses 
of both DFT and of model Hamiltonian approaches. The 
figure shows the exact ground state density of ten atoms 
with interatomic spacing 6 = 4. The edges induce a 
staggered pattern of strong and weak bonds which de- 
cays slowly into the bulk, and is therefore significant 
throughout this small system. Qualitatively, we can un- 
derstand the staggered behavior from a 10 site Hubbard 
model (at half filling with U/t — 4 chosen arbitrarily) 
or a 10 site Heisenberg model. The Heisenberg ground 
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FIG. 3. Entanglement spectrum at the center of interacting 
4-atom chains with various interatomic separations b. Nl 
refers to the number of electrons to the left of the cut for 
each density matrix eigenstate. 
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FIG. 2. Spontaneous dimerization of the density for a chain 
of 10 soft hydrogen atoms with interatomic spacing b — 4 
(dashed lines are a guide to the eye). The upper panel com- 
pares the densities predicted by DFT within the LDA ap- 
proximation; the lower panel shows the spin densities for the 
unrestricted LDA. Also shown is the expectation value of the 
kinetic energy for a Hubbard model and the exchange energy 
for the Heisenberg model on 10 lattice sites. 



state has resonating valence bond character; in a per- 
fect near-neighbor RVB state, the edges would suppress 
all resonance and drive the weak bonds to zero. The ac- 
tual Heisenberg ground state has longer range resonances 
which reduce these effects. In the Hubbard model, the 
strong exchange bonds show up as bonds with lower ki- 
netic energy. However, neither lattice model reveals the 
increased electron density on the strong bonds, stemming 
from the strong hopping. The exact results suggest that 
these models might be improved by bond-dependent in- 
teractions t and J. The DFT calculations capture even 
fewer properties of the true ground state. Unrestricted 
LDA predicts an energy £\jlda = —11.364 which is close 
to the exact energy Eo = —11.496, but it predicts no 
staggered bond density and breaks spin symmetry, pro- 
ducing a long-ranged antiferromagnetic state as shown 
in the lower panel of Fig. 2. Restricted LDA captures 
the staggered density pattern qualitatively, but gives a 
slightly higher energy -Erlda = —11.323 and fails to re- 
produce the correct local spin correlations since its wave- 
function is a Slater determinant of extended orbitals. 

The onset of strong correlation with increasing bond 
length is often identified with the Coulson-Fischer point 
[14], where an unrestricted Hartree-Fock calculation 
spontaneously breaks spin symmetry. A different way to 
distinguish strong from weak correlation is by examining 
the entanglement spectrum, readily accessible in DMRG. 
Defining the left reduced density matrix pl — Ttr|\I/) (\I/|, 
where the trace is over all grid points in the right half of 
the system, the entanglement spectrum consists of the 
energies of the entanglement Hamiltonian He = — In pi, 
[15]. The most probable density matrix eigenstates are 
those in the low 'energy' part of the entanglement spec- 
trum. By classifying these states according to their total 



particle number Nl, we can understand the dominant 
quantum fluctuations of the ground state. Figure 3 shows 
the entanglement spectrum at the center of a series of 
four-atom chains with increasing interatomic separation. 
A sharp crossover at b ~ 6, where the probability for 
charge fluctuations drops below that of pure spin fluctu- 
ations, signals the onset of strongly correlated behavior. 
The foundation of DFT is the Hohenberg-Kohn theo- 
rem, which states that for a specific type of interaction 
and a ground state electron density n(x), the potential 
v(x) giving this density is unique if it exists [16]. This 
implies that the total energy and various terms such as 
the kinetic energy are functionals of the density. Stan- 
dard DFT approaches are based on approximate func- 
tionals where the approximations are explicit formulas. 
However, having an exact many-body Schrodinger equa- 
tion solver make it possible to implement the functionals 
exactly, not as formulas, but as algorithms. The key 
step in each algorithm is performing an inversion, an it- 
erative calculation of the potential which gives rise to 
a specified density. (Note that an interacting inversion 
is much more difficult than the non-interacting inversion 
needed to compute the Kohn-Sham potential from the ex- 
act ground state density [17, 18].) Interacting inversions 
have previously been done only for a Id two-electron sys- 
tem during a time evolution [12]. To demonstrate that 
we can perform interacting inversions of any reasonable 
density (we exclude densities that vary appreciably on 
the grid scale, for instance), Fig. 4 shows an interacting 
inversion of the density h(x) = e x / 1 5-x-/2+x /w-x /750 
normalized to contain four electrons. We begin with a 



trial potential v(x) from which we obtain the trial den- 
sity n(x) using DMRG. Then the potential is updated 
as v(x) — > v(x) + a(n{x) — n(x)) for some a > 0. This 
procedure can be motivated heuristically or from a min- 
imization principle. By taking a — 0.5, the inversion in 
Fig. 4 converged in 300 steps such that the maximum 
difference between the inverted and target densities was 
less than 0.01. Once an inversion is done, it is simple to 
obtain exact density functionals, such as the ground state 
energy ^[n], because one has obtained the many-body 
wavefunction. 

Many oxide materials of current interest are too 
strongly correlated for present DFT methods, but cru- 
cial properties must be calculated to an accuracy far be- 
yond that of simple model Hamiltonians. The method de- 
scribed here provides a new, alternative route to studying 
strongly correlated systems. All existing approximations, 
from heuristic corrections to standard functionals, such 
as LDA+U [19], to methods developed for lattice models, 
such as dynamical mean field theory [20] , can be applied 
and tested more easily, thoroughly, and accurately in the 
present setting. Because our Id world captures a feature 
crucial to density functional approximations, namely the 
continuum instead of a lattice, such studies should pro- 
vide the insight needed to construct more accurate den- 
sity functionals for real strongly-correlated materials. 
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FIG. 4. Interacting inversion of an artificial density for a four 
electron system, labeled here as the target density. The lower 
plot shows trial potentials at 5 selected steps. The dashed 
curves in the upper plot are the corresponding trial densities 
found by computing the ground state in each potential with 
interactions given by v cc (x) = 1/y/x 2 + 1. 
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